Within this chapter developmental plasticity was explored within Acanthochromis polyacanthus that were collected from two different regions (i.e., low-latitude, Cairns, and high-latitude, Mackay). Fish were held in common garden experiments at 28.5 C. Clutch data was collected along with parental morphometeric data to determine how fish from each population performed at 28.5 C, a temperature that was shown to produce similar absoluate aerobic scope performances in a previous study [link]. Once hatched offspring were placed into three different temperature treatments, 28.5, 30, and 31.5 C. After approximately 5-6 months offspring length and weight was measured, as well as CTmax and respiration during CTmax trials.
library(tidyverse) # data manipulation
library(ggpubr) # figure arrangement
library(brms) # Bayesian models
library(StanHeaders)# needed to run Bayesian models
library(rstan) # needed to run Bayesian models
library(standist) # needs to be installed
library(bayesplot) # needed for MCMC diagnostics
library(DHARMa) # model validation
library(ggdist) # partial plots
library(tidybayes) # partial plots
library(broom.mixed) # model investigation
library(emmeans) # pairwise comparisons
library(rstanarm) # pairwise comparisons - need for emmeans
library(ggpubr) # arranging figures
library(mgcv) # GAM
library(modelr) # plotting
library(segmented) # piece wise regression
library(strucchange) #piece wise regressionThis data set has no point outliers however, in the next chunk 3 samples will be discarded due to poor data quailty
rows_of_data <- count(lat_resp_dat, sampleID)
lat_resp_dat2 <- lat_resp_dat |>
mutate(dev.temp = as.factor(dev.temp),
replicate = str_sub(sampleID, -1,-1),
population = factor(population)) |>
# number of observations = 5758
filter(sampleID != "56.CARL.137.28,5.1", # 5777 - 76 = 5682
sampleID != "56.CARL.137.28,5.2", # 5701 - 64 = 5618
sampleID != "60.LCKM.152.30.1" # 5637 - 76 = 5542
) |>
filter(time_lag_sec >2001) |> # remove samples from first 5 cycles (i.e., first 33 minutes)
group_by(sampleID) |>
mutate(max_value_index = which.max(rate.output),
row_number = row_number(),
max.rate = max(rate.output),
low.rate = mean(rate.output[order(rate.output)[1:10]]),
net.rate = max.rate - low.rate) |>
filter(row_number <= max_value_index) |>
ungroup() |>
mutate(region = factor(region),
dev.temp =factor(dev.temp),
level = as.factor(case_when(tank >= 1 & tank <= 199 ~ 1,
tank >= 200 & tank <= 299 ~ 2,
tank >= 300 & tank <= 399 ~ 3,
TRUE ~ NA_real_)),
female = factor(female)) |>
unite("dev.temp.region",c(dev.temp,region), remove=FALSE) |>
mutate(dev.temp.region = factor(dev.temp.region))eda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 2, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.figeda1 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output)) +
geom_point(alpha=0.5, color = "black") +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
theme_classic() +
ggtitle("All data points - 2nd order polynomial")
eda2 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda3 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ region) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda4 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=region)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ dev.temp) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Region - 2nd order polynomial")
eda5 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=population)) +
geom_point(alpha=0.5) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE), color = "red") +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda6 <- ggplot(lat_resp_dat2, aes(x=Temperature, y=rate.output, color=dev.temp)) +
geom_point(alpha=0.1) +
geom_smooth(method="lm",
se=TRUE,
fill=NA,
formula = y~poly(x, 3, raw=TRUE)) +
facet_wrap(~ population) + theme_classic() +
theme(legend.position = "none") +
ggtitle("Population - 2nd order polynomial")
eda.fig <- ggarrange(eda1,eda2,eda3,eda4,eda5,eda6,
nrow = 6,
ncol=1); eda.figFirst we will place random factors only within out model and see if they do a good job explaining the variance within our model. We will also be include priors which will be based off of out length data.
Hypothesis test will include:
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
f.model.null <- bf(rate.output ~ 1,
family = gaussian())
model.null <- brm(f.model.null,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
f.model1 <- bf(rate.output ~ 1 + (1| female) + (1| tank),
family=gaussian())
model1 <- brm(f.model1,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
f.model2 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level),
family=gaussian())
model2 <- brm(f.model2,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 13 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
f.model3 <- bf(rate.output ~ 1 + (1| female) + (1| tank) + (1| level) + (1| population),
family=gaussian())
model3 <- brm(f.model3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 5 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
f.model4 <- bf(rate.output ~ 1 + (1| female) + (1 | tank) + (1| level) + (1| population)+ (1| clutch_order),
family=gaussian())
model4 <- brm(f.model4,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Output of model 'model.null':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -31469.1 78.1
## p_loo 2.5 0.4
## looic 62938.2 156.2
## ------
## MCSE of elpd_loo is 0.0.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.0]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.6 99.6
## p_loo 57.6 3.1
## looic 61325.1 199.1
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.7 99.5
## p_loo 57.6 3.0
## looic 61325.4 198.9
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.7 99.6
## p_loo 58.1 3.2
## looic 61325.5 199.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model4':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30662.5 99.4
## p_loo 57.5 3.0
## looic 61325.0 198.8
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model4 0.0 0.0
## model1 -0.1 0.4
## model2 -0.2 0.4
## model3 -0.2 0.4
## model.null -806.6 48.2
lat_resp_dat2 |>
group_by(region,dev.temp) |>
summarise(mean(na.omit(rate.output)),
sd(na.omit(rate.output))) ## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.
rate.priors <- prior(normal(430, 0.25), class="Intercept") +
prior(normal(0, 40), class="b")
prior(student_t(3, 0, 195), class = "sigma")f.model1.p1 <- bf(rate.output ~ 1 +
dev.temp*region + Temperature +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p1 <- brm(f.model1.p1,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p1, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p1, ndraws=250, summary=FALSE)
model1.p1_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p1_resids) ; testDispersion(model1.p1_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98596, p-value = 0.704
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
f.model1.p2 <- bf(rate.output ~ 1 +
dev.temp*region + poly(Temperature, 2) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p2 <- brm(f.model1.p2,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p2, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p2, ndraws=250, summary=FALSE)
model1.p2_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p2_resids) ; testDispersion(model1.p2_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98971, p-value = 0.688
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
f.model1.p3 <- bf(rate.output ~ 1 +
dev.temp*region + poly(Temperature, 3) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank) ,
family=gaussian())
model1.p3 <- brm(f.model1.p3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
#--- distribution check ---#
pp_check(model1.p3, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.p3, ndraws=250, summary=FALSE)
model1.p3_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.p3_resids) ; testDispersion(model1.p3_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.98512, p-value = 0.616
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
## Output of model 'model1.p1':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -29946.0 94.1
## p_loo 60.3 3.1
## looic 59892.1 188.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.1]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p2':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30525.1 100.7
## p_loo 58.9 3.2
## looic 61050.1 201.3
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.7, 1.2]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Output of model 'model1.p3':
##
## Computed from 1800 by 4776 log-likelihood matrix.
##
## Estimate SE
## elpd_loo -30516.1 100.1
## p_loo 58.7 3.1
## looic 61032.1 200.2
## ------
## MCSE of elpd_loo is 0.2.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.8, 1.2]).
##
## All Pareto k estimates are good (k < 0.69).
## See help('pareto-k-diagnostic') for details.
##
## Model comparisons:
## elpd_diff se_diff
## model1.p1 0.0 0.0
## model1.p3 -570.0 31.3
## model1.p2 -579.0 31.5
f.model1.gamk4 <- bf(rate.output ~ 1 +
t2(dev.temp,region,Temperature, k=4, bs="fs", by=dev.temp.region) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
model1.gamk4 <- brm(f.model1.gamk4,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 1 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
#--- distribution check ---#
pp_check(model1.gamk4, type = 'dens_overlay_grouped', ndraws=150, group="dev.temp.region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk4, ndraws=250, summary=FALSE)
model1.gamk4_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.gamk4_resids) ; testDispersion(model1.gamk4_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97661, p-value = 0.472
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
### {-}
f.model1.gamk3 <- bf(rate.output ~ 1 +
t2(dev.temp,region,Temperature, k=3, bs="fs", by=dev.temp.region) +
scale(mass, center=TRUE, scale=TRUE) +
(1| female) + (1| tank),
family=gaussian())
model1.gamk3 <- brm(f.model1.gamk3,
data = lat_resp_dat2,
prior = rate.priors,
warmup = 500,
iter = 5000,
seed=123,
cores=2,
save_pars = save_pars(all=TRUE),
sample_prior = "yes",
chains = 2,
thin = 5,
control = list(adapt_delta=0.95)) ## Compiling Stan program...
## Start sampling
## Warning: There were 3 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
#--- distribution check ---#
pp_check(model1.gamk3, type = 'dens_overlay_grouped', ndraws=150, group="region")#--- DHARMa checks ---#
preds <- posterior_predict(model1.gamk3, ndraws=250, summary=FALSE)
model1.gamk3_resids <- createDHARMa(simulatedResponse = t(preds),
observedResponse = lat_resp_dat2$rate.output,
fittedPredictedResponse = apply(preds, 2, mean),
integerResponse = 'student')
plot(model1.gamk3_resids) ; testDispersion(model1.gamk3_resids)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.97942, p-value = 0.448
## alternative hypothesis: two.sided
## 'pars' not specified. Showing first 10 parameters by default.
## 'pars' not specified. Showing first 10 parameters by default.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## 'pars' not specified. Showing first 10 parameters by default.
### {-}
## Warning in AIC.default(model1.gamk3, model1.gamk4): models are not all fitted
## to the same number of observations
## Warning: There were 3 divergent transitions after warmup. Increasing
## adapt_delta above 0.95 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: rate.output ~ 1 + t2(dev.temp, region, Temperature, k = 4, bs = "fs", by = dev.temp.region) + scale(mass, center = TRUE, scale = TRUE) + (1 | female) + (1 | tank)
## Data: lat_resp_dat2 (Number of observations: 4776)
## Draws: 2 chains, each with iter = 5000; warmup = 500; thin = 5;
## total post-warmup draws = 1800
##
## Smooth Terms:
## Estimate Est.Error
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 420.23 192.27
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 409.00 203.52
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 807.91 370.10
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 838.69 338.93
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 672.50 274.67
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 685.97 275.46
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 910.31 389.52
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 989.69 425.40
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 423.24 186.68
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 454.59 194.73
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 592.65 257.89
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 543.84 236.79
## l-95% CI u-95% CI
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 204.95 935.83
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 166.01 941.29
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 392.88 1758.28
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 413.66 1700.17
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 350.27 1420.75
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 338.90 1390.80
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 457.39 1876.23
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 482.79 2114.40
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 202.85 861.71
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 217.51 970.44
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 291.71 1249.28
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 264.18 1129.12
## Rhat Bulk_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 1.00 1819
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 1.00 1648
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 1.00 1630
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 1.00 1600
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 1.00 1741
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 1.00 1751
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 1.00 1776
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 1.00 1760
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 1.00 1883
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 1.00 1700
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 1.00 1599
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 1.00 1583
## Tail_ESS
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_1) 1829
## sds(t2dev.tempregionTemperaturedev.temp.region28_core_2) 1746
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_1) 1529
## sds(t2dev.tempregionTemperaturedev.temp.region28_leading_2) 1662
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_1) 1761
## sds(t2dev.tempregionTemperaturedev.temp.region30_core_2) 1706
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_1) 1714
## sds(t2dev.tempregionTemperaturedev.temp.region30_leading_2) 1744
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_1) 1911
## sds(t2dev.tempregionTemperaturedev.temp.region31_core_2) 1708
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_1) 1606
## sds(t2dev.tempregionTemperaturedev.temp.region31_leading_2) 1563
##
## Group-Level Effects:
## ~female (Number of levels: 15)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 118.65 31.65 65.99 189.82 1.00 1128 1388
##
## ~tank (Number of levels: 52)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 82.36 10.87 64.11 107.11 1.00 1276 1524
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat
## Intercept 430.00 0.24 429.54 430.47 1.00
## scalemasscenterEQTRUEscaleEQTRUE 19.55 3.77 12.02 26.65 1.00
## Bulk_ESS Tail_ESS
## Intercept 1801 1707
## scalemasscenterEQTRUEscaleEQTRUE 1683 1613
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 114.35 1.19 111.99 116.71 1.00 1771 1784
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#model1.re.wo |> get_variables()
model1.gamk4 |> gather_draws(`b_.*|sigma`, regex =TRUE) |>
median_hdci()model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="dev.temp") |> summary(infer=TRUE)## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)
model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)
## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)
## region dev.temp emmean lower.HPD upper.HPD
## core 28 444 367 511
## leading 28 419 303 511
## core 30 467 399 538
## leading 30 399 295 512
## core 31 477 400 547
## leading 31 485 393 577
##
## Point estimate displayed: median
## HPD interval probability: 0.95
model1.gamk4 |> emmeans(pairwise ~ region*dev.temp, type="response") |> pairs(by="region") |> confint()## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)
## NOTE: A nesting structure was detected in the fitted model:
## dev.temp.region %in% (dev.temp*region)
## $emtrends
## dev.temp region Temperature.trend lower.HPD upper.HPD
## 28 core 44.872 31.78 57.1
## 30 core 2.862 -8.80 15.7
## 31 core 24.090 13.29 35.6
## 28 leading -0.815 -14.74 12.8
## 30 leading 1.168 -13.99 16.4
## 31 leading 7.000 -5.06 19.1
##
## Point estimate displayed: median
## HPD interval probability: 0.95
##
## $contrasts
## contrast estimate lower.HPD upper.HPD
## dev.temp28 core - dev.temp30 core 42.04 24.173 60.14
## dev.temp28 core - dev.temp31 core 20.95 4.017 38.19
## dev.temp28 core - dev.temp28 leading 45.62 26.541 63.58
## dev.temp28 core - dev.temp30 leading 43.94 25.948 64.08
## dev.temp28 core - dev.temp31 leading 38.00 20.584 55.40
## dev.temp30 core - dev.temp31 core -21.32 -37.216 -3.77
## dev.temp30 core - dev.temp28 leading 3.65 -14.587 23.10
## dev.temp30 core - dev.temp30 leading 1.71 -16.474 22.04
## dev.temp30 core - dev.temp31 leading -4.18 -20.285 13.35
## dev.temp31 core - dev.temp28 leading 24.81 6.454 42.12
## dev.temp31 core - dev.temp30 leading 23.23 5.496 42.63
## dev.temp31 core - dev.temp31 leading 17.00 -0.205 32.99
## dev.temp28 leading - dev.temp30 leading -1.55 -20.579 19.27
## dev.temp28 leading - dev.temp31 leading -7.62 -25.674 10.96
## dev.temp30 leading - dev.temp31 leading -6.18 -25.762 12.80
##
## Point estimate displayed: median
## HPD interval probability: 0.95
resp.plot <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~dev.temp); resp.plotresp.plot2 <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic() +
facet_wrap(~region); resp.plot2resp.plot3 <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_fitted_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
ggplot(aes(x=Temperature, color=exp_group)) +
geom_line(aes(y=.value, group=paste(exp_group, .draw)), alpha=0.05) +
geom_smooth(aes(y=.value, color=exp_group), se=FALSE, size=1.5) +
theme_classic(); resp.plot3posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "28_core")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 70
## m = 2 60 81
## m = 3 54 70 85
## m = 4 40 55 70 85
## m = 5 15 40 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.7
## m = 2 0.6 0.81
## m = 3 0.54 0.7 0.85
## m = 4 0.4 0.55 0.7 0.85
## m = 5 0.15 0.4 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 1170431.8 216925.0 88460.8 49097.0 44088.4 41482.4
## BIC 1229.8 1070.4 989.9 940.3 938.7 941.8
## [1] 33.86963 35.04506 36.14703
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.86963 35.04506 36.14703
my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.86963 32.14501 0.02701804
## psi2.x 35.04506 33.52461 0.02942867
## psi3.x 36.14703 34.93830 0.03209560
## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -15.086 0.53064 -28.429 -16.139 -14.032
## slope2 25.029 1.05370 23.754 22.937 27.122
## slope3 65.828 1.05370 62.474 63.736 67.921
## slope4 99.259 0.48165 206.080 98.303 100.220
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "30_core")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 78
## m = 2 70 85
## m = 3 55 70 85
## m = 4 15 55 70 85
## m = 5 15 30 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.78
## m = 2 0.7 0.85
## m = 3 0.55 0.7 0.85
## m = 4 0.15 0.55 0.7 0.85
## m = 5 0.15 0.3 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 1289086.6 223458.9 95112.0 79870.4 72776.6 72689.6
## BIC 1239.4 1073.4 997.2 988.9 988.8 997.9
## [1] 33.86963 35.04506 36.14703
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.86963 35.04506 36.14703
my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.86963 31.47574 0.05887590
## psi2.x 35.04506 34.20755 0.02233563
## psi3.x 36.14703 35.37701 0.02265965
## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 26.77600 1.37500 19.4730 24.04500 29.507
## slope2 0.81219 0.58750 1.3825 -0.35463 1.979
## slope3 75.59800 2.06930 36.5340 71.48800 79.708
## slope4 154.35000 0.99772 154.7000 152.37000 156.330
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic() ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "31_core")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 74
## m = 2 64 83
## m = 3 54 70 85
## m = 4 40 55 70 85
## m = 5 22 40 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.74
## m = 2 0.64 0.83
## m = 3 0.54 0.7 0.85
## m = 4 0.4 0.55 0.7 0.85
## m = 5 0.22 0.4 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 951373.7 176567.4 73635.2 44061.0 39736.9 39265.8
## BIC 1209.0 1049.8 971.6 929.4 928.3 936.3
## [1] 33.94309 35.04506 36.14703
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 33.16127 0.02769747
## psi2.x 35.04506 34.42367 0.02225931
## psi3.x 36.14703 35.50138 0.02594152
## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 5.4424 0.21502 25.311 5.0153 5.8694
## slope2 33.1200 0.89665 36.937 31.3390 34.9010
## slope3 75.5130 1.08240 69.767 73.3630 77.6630
## slope4 108.9600 0.53408 204.020 107.9000 110.0300
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "28_leading")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 79
## m = 2 70 85
## m = 3 55 70 85
## m = 4 15 55 70 85
## m = 5 15 38 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.79
## m = 2 0.7 0.85
## m = 3 0.55 0.7 0.85
## m = 4 0.15 0.55 0.7 0.85
## m = 5 0.15 0.38 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 1612756 264488 108644 99987 98860 97905
## BIC 1262 1090 1010 1011 1019 1028
## [1] 33.94309 35.04506 36.14703
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 34.20914 0.05958793
## psi2.x 35.04506 34.93752 0.07430573
## psi3.x 36.14703 35.75777 0.08572172
## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 1.386 0.61625 2.2490 0.16205 2.6099
## slope2 66.933 8.64990 7.7379 49.75300 84.1120
## slope3 135.080 7.49110 18.0320 120.20000 149.9600
## slope4 186.000 2.83140 65.6920 180.37000 191.6200
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "30_leading")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 78
## m = 2 70 85
## m = 3 55 70 85
## m = 4 15 55 70 85
## m = 5 15 37 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.78
## m = 2 0.7 0.85
## m = 3 0.55 0.7 0.85
## m = 4 0.15 0.55 0.7 0.85
## m = 5 0.15 0.37 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 2321698 379911 154479 140940 140332 138992
## BIC 1298 1126 1046 1046 1054 1063
## [1] 33.94309 35.04506 36.14703
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 34.13619 0.04916754
## psi2.x 35.04506 34.92604 0.06043609
## psi3.x 36.14703 35.75675 0.07509539
## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 0.24934 0.64634 0.38578 -1.0343 1.533
## slope2 78.34700 7.65450 10.23500 63.1450 93.550
## slope3 161.85000 7.65450 21.14400 146.6400 177.050
## slope4 221.18000 2.89310 76.44900 215.4300 226.920
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
posterior_draws <- lat_resp_dat2 |>
group_by(region, dev.temp,dev.temp.region) |>
data_grid(Temperature=seq(from=min(lat_resp_dat2$Temperature),
to=max(lat_resp_dat2$Temperature),
length.out=100),
mass=0.00148299) |>
add_linpred_draws(model1.gamk4,
n=300,
re_formula=NA) |>
unite("exp_group", c("region","dev.temp"), remove=FALSE) |>
filter(dev.temp.region == "31_leading")
spline_fit <- smooth.spline(posterior_draws$Temperature, posterior_draws$.linpred, nknots = 4)
ggplot(posterior_draws, aes(x=Temperature, y=.linpred)) +
geom_point(color="black", alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y,
color="red"),
size=2) ## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
spline.data <- as.data.frame(x= spline_fit$x) |>
mutate(y = spline_fit$y)|>
rename(x = `spline_fit$x`)
plot(spline.data)##
## Optimal (m+1)-segment partition:
##
## Call:
## breakpoints.formula(formula = spline.data$y ~ 1)
##
## Breakpoints at observation number:
##
## m = 1 78
## m = 2 70 85
## m = 3 55 70 85
## m = 4 26 55 70 85
## m = 5 17 32 55 70 85
##
## Corresponding to breakdates:
##
## m = 1 0.78
## m = 2 0.7 0.85
## m = 3 0.55 0.7 0.85
## m = 4 0.26 0.55 0.7 0.85
## m = 5 0.17 0.32 0.55 0.7 0.85
##
## Fit:
##
## m 0 1 2 3 4 5
## RSS 995002.7 164741.2 67258.5 60638.9 57568.3 57272.2
## BIC 1213.5 1042.9 962.5 961.4 965.4 974.1
## [1] 33.94309 35.04506 36.14703
Use out from strucchange summary to number of breakpoints as well as approximate location of breakpoints
## [1] 33.94309 35.04506 36.14703
my.seg <- segmented(my.lm,
seg.Z = ~x,
psi=list(x=spline.data[breakpoint.values, c("x")]))
#summary(my.seg)
my.seg$psi ; slope(my.seg)## Initial Est. St.Err
## psi1.x 33.94309 33.54958 0.02768565
## psi2.x 35.04506 34.58443 0.02256263
## psi3.x 36.14703 35.58101 0.02554307
## $x
## Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -7.8319 0.24673 -31.743 -8.3219 -7.3418
## slope2 31.0790 1.61930 19.193 27.8630 34.2960
## slope3 91.3370 1.61930 56.403 88.1200 94.5530
## slope4 138.6300 0.76778 180.560 137.1100 140.1600
my.fitted <- fitted(my.seg)
my.model <- data.frame(temperature = spline.data$x, rate.output= my.fitted)
ggplot(my.model, aes(x=temperature, y=rate.output)) + geom_line()ggplot(posterior_draws, aes(x=Temperature, y=.linpred, color=.linpred)) +
geom_point(alpha=0.05, size=1) +
geom_smooth(data=as.data.frame(spline_fit$data), aes(x=x, y=y), color="red",
size=3) + ylab("rate.output")+
geom_vline(xintercept = my.seg$psi[1,2]) +
geom_vline(xintercept = my.seg$psi[2,2])+
geom_vline(xintercept = my.seg$psi[3,2])+
geom_vline(xintercept = my.seg$psi[1,1], color="blue") +
geom_vline(xintercept = my.seg$psi[2,1], color="blue") +
geom_vline(xintercept = my.seg$psi[3,1], color="blue") +
geom_line(data=my.model, aes(x=temperature, y=rate.output), color="black", size=2) +
theme_classic()## `geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'